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ABSTRACT 



Context. We study the /c-mechanism that excites radial oscillations in Cepheid variables. 

Aims. We address the mode couplings that manage the nonlinear saturation of the instability in direct numerical simulations (DNS). 

Methods. We project the DNS fields onto an acoustic subspace built from the regular and adjoint eigenvectors that are solutions to 

the linear oscillation equations. 

Results. We determine the time evolution of both the amplitude and kinetic energy of each mode that propagates in the DNS. More 

than 98% of the total kinetic energy is contained in two modes that correspond to the linearly-unstable fundamental mode and the 

linearly- stable second overtone. Because the eigenperiod ratio is close to 1/2, we discover that the nonlinear saturation is due to a 2:1 

resonance between these two modes. An interesting application of this result concerns the reproduction of Hertzsprung 's progression 

observed in Bump Cepheids. 
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1. Introduction 

In Gastine & Dintrans (2008) (hereafter Paper I), we modelled 
the A:-mechanism in classical Cepheids using a simplified ap- 
proach, that is, the propagation of acoustic waves in a layer of 
gas where the opacity bump is represented by a hollow in ra- 
diative conductivity. To maintain the mode excitation, the two 
following conditions apply: 

- A sufficient width and amplitude are required for the con- 
ductivity hollow. 

- This hollow must be located in a precise zone called the 
"transition region". 

Our parametric approach enabled us to check the reality of 
these two conditions and to obtain the instability strips from a 
linear- stability analysis. By starting from the most favourable 
situations (i.e. the most unstable linear modes), we performed 
the corresponding 1-D and 2-D nonlinear Direct Numerical 
Simulations (hereafter DNS). These DNS confirmed with a note- 
worthy agreement both the growth rates and spatial structures of 
linearly-unstable modes. The nonlinear saturation was reached 
and a quantitative study of the involved mode couplings re- 
mained to be done. 

This is the principle aim of this paper. To study the nonlin- 
ear interactions between modes, we adapt a powerful method 
used in, e.g., aeroacoustics (Salwen & Grosch 1981; Wang et al. 
2006) or numerical astrophysics (Bogdan et al. 1993, hereafter 
BCM93). It is based on a projection of DNS fields onto a basis 
shaped from the regular and adjoint eigenvectors that are solu- 
tions to the linear oscillation equations. This method has already 
been used in a simplified version by one of us to study internal 
waves in DNS (Dintrans & Brandenburg 2004; Dintrans et al. 
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2005). In the quoted investigation the eigenvalue problem was 
adiabatic, and therefore Hermitian, and the authors only ac- 
counted for projections onto regular eigenf unctions. This reduc- 
tion is not applicable in our A:-mechanism simulations since the 
transition region requires low densities close to the surface and 
then high difTusivities (Paper I). Eigenmodes are highly non- 
adiabatic and imply non-Hermitian oscillation operators with 
non-orthogonal eigenf unctions. Solving the adjoint problem is 
therefore mandatory to determine both mode amplitudes and en- 
ergy. 

By using projections onto the two respective sets of eigen- 
vectors, the regular and adjoint ones, the time evolution of each 
acoustic mode propagating in DNS is completed. The kinetic en- 
ergy in each acoustic mode is available that highlights the energy 
transfer between modes. One of the main results of this work is 
that more than 98% of the total kinetic energy is contained in 
both the fundamental mode and the second overtone. Because 
this second overtone is linearly stable, we show that its large am- 
plitude results from a 2: 1 resonance with the fundamental mode. 
In our numerical experiments, this resonance is efficient because 
the ratio of the two involved periods, that is P2/P0, is close to 
1/2. 

This nonlinear saturation based on a 2:1 resonance is an in- 
teresting result because this mechanism has been proposed to 
explain the secondary bump in light curves of some classical 
Cepheids named "Bump Cepheids". The bump position is cor- 
related with the oscillation period that leads to the well-known 
Hertzsprung progression (Hertzsprung 1926; Payne- Gaposchkin 
1954). The bump first appears on the descending branch of the 
light curve of Population I Cepheids with periods of about 6- 
7 days and then travels up this curve to reach its maximum for 
10-11 day Cepheids. For longer periods, it moves down in the as- 
cending branch and disappears for periods longer than 20 days 
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(Whitney 1983; Tsvetkov 1990; Bono et al. 2000). Two main 
theories have been developed to explain this phenomenon: 

- Whitney (1956) suggested an "echo mechanism" that 
Christy (1968) developed in more detail. In this model, two 
radial-pressure waves (called Christy's waves) moving in op- 
posite directions are generated during compression phases in 
the Hell ionisation region. The pulse running downward re- 
flects on the stellar core and reaches the surface at the next 
period during which the bump appears. 

- The second explanation, known as the "resonance mech- 
anism", has been proposed by Simon & Schmidt (1976) 
(hereafter SS76). The bump results from a resonance be- 
tween the fundamental mode and the second overtone that 
is possible for a period ratio P2/P0 close to 1/2. 

Three papers of Whitney (1983) and Aikawa & Whitney (1984- 
1985) compared the echo and resonance theories. The basic 
idea was to consider both approaches as two complementary 
sides of the same physical process. Despite promising results 
on polytropes (Whitney 1983), the dynamical approach showed 
that the acoustic-ray formalism did not reproduce the "run- 
ning waves" in Christy's diagram. This implies that most of the 
energy is contained in standing waves rather than progressive 
waves (Aikawa & Whitney 1984, 1985). 

The resonance mechanism has also been investigated us- 
ing the amplitude-equations formalism by Buchler & Goupil 
(1984) and Klapp et al. (1985). They established the domi- 
nant role of a 2:1 resonance because phases and amplitudes 
of modes show characteristic variations with the period-ratio 
P2/P0 (see also Simon & Lee 1981; Kovacs & Buchler 1989; 
Buchler et al. 1990). Despite these results favoring the resonance 
phenomenon, the echo-resonance controversy remains topical 
(Fadeyev & Muthsam 1992; Bono et al. 2000, 2002). 

It is interesting to know whether our model is able to re- 
produce this Hertzsprung progression. Up to 400 DNS are com- 
pleted to cover a significant range in the P2/P0 ratio, while 
studying the bump position in luminosity curves. This allows 
us to accurately reproduce the expected bump-progression with 
the ratio value: (i) if P2/P0 > 1/2, the bump is located in the 
descending branch; (ii) if P2/P0 < 1/2, the bump is located in 
the ascending branch. 

In Sects. 2 and 3, we introduce the general-oscillations equa- 
tions and the associated adjoint problem, respectively. We de- 
velop our projection method and provide results in Sect. 4. 
Sect. 5 concerns the Hertzsprung progression and we outline our 
conclusions in Sect. 6. 



2. The pulsation model and corresponding DNS 

We focus on radial modes propagating in Cepheids and re- 
strict our study to the 1-D case. Our system, which represents 
a local zoom around an ionisation region, is composed of a 
layer of width d filled with a monatomic and perfect gas with 
y = c p lc v = 5/3 (c p and c v are specific heats). Gravity g = -ge z 
and kinematic viscosity v are assumed to be constant. Following 
Paper I, the ionisation region is described by a parametric hollow 
in radiative conductivity that corresponds to a bump in opac- 
ity. We recall below the temperature-dependent profile that is 
adopted for the radiative conductivity 
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where 7b U mp is the position of the hollow in temperature and <x, 
e and ^R denote its slope, width and relative amplitude, respec- 
tively. 

2.1. Instability strips from the linear-stability analysis 

We are interested in small perturbations about the hydrostatic 
and radiative equilibria. The layer is fully radiative and the dif- 
fusion approximation implies that 



f' = -K vr - K'VT 



(3) 



for the radiative flux perturbation (hereafter the "0" subscripts 
mean equilibrium quantities, while primes denote Eulerian 
ones). 

Following Paper I, the depth d of the layer is selected as the 
length scale, and the top density p top and temperature 7\ op as den- 
sity and temperature scales, respectively. The velocity scale is 
thus (c p T i0V ) 112 . Finally, gravity and fluxes are provided in units 
of c p T top /d andp t0 p(c /? r top ) 3/2 , respectively. The linearised per- 
turbations obey the following temperature, momentum and con- 
tinuity dimensionless equations: 
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AR = u, 

dz dz 

where R = p'/po denotes the density perturbation, u the velocity 
and fbot the imposed bottom flux. Here we seek normal modes 
using a time-dependence of the form exp(At), with A = r + id) (r 
is the damping or growth rate of the mode and co its frequency). 
Finally, by assuming rigid walls at both limits of the domain 
for the velocity, a perfect conductor at the bottom and a perfect 
insulator at the top for the temperature, we obtain the following 
boundary conditions 

M = 0forz = (0,l), 



dV 



= for z = 0, 



dz 
[T f = 0forz= 1. 
Eqs. (4-5) define an eigenvalue problem 

Ax w = A n x n , 



(5) 



(6) 



where A n is the eigenvalue associated with the (complex) eigen- 
vector *P„ = (7", u, R) T . Here n defines the mode order, that is, 
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e = 0.4 and a = 7 




n 



- o 

i 



2.8 2.6 2.4 2.2 2.0 1.8 1.6 



Table 1. First linear eigenvalues of the unstable setup used to 
start the DNS (this setup is denoted by the white cross in the 
instability strip in Fig. 1). Note that only the fundamental n = 
mode is linearly unstable. 
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0.361 
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-9.363 


0.304 
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-11.91 


0.263 



^bump 



Fig. 1. Instability strip for the fundamental mode in the plane 
(Tbump, &) for e = 0.4, cr = 7 and v = 1 x 10 -4 . The top (at 
z = 1) density and temperature are p = 2.5 x 10" 3 and T = 1, 
respectively, while gravity is g = 7 everywhere. The white cross 
denotes a particular simulation further studied for 7b um p = 2.1 
andJ?l = 70%. 



the number of nodes of the eigenfunction u while the matrix A is 
a differential operator provided in Eq. (7), where the differential 
notation D = d/dz is used for clarity. 

Figure 1 displays the instability strip for the fundamen- 
tal mode 1 . In this parametric survey, we fix the slope cr and 
width e of the conductivity hollow, whereas its amplitude J{ 
and temperature 7b U mp vary (Eq. 1). For every couple (7b um p, &), 
both equilibrium fields and solutions to the eigenvalue problem 
(6) are completed using the LSB code (Linear Solver Builder, 
Valdettaro et al. 2007). Unstable modes with r > are identified 
among all eigenvalues in each spectrum and only the fundamen- 
tal mode is excited by the A:-mechanism in these simulations. 

2.2. The associated nonlinear DNS 

To confirm the reality of the instability strip and study the mode 
saturation, we perform a DNS of the nonlinear problem. We start 
from a linearly-unstable setup discovered in the parametric sur- 
vey (the white cross in Fig. 1 of which the oscillation spectrum 
is provided in Table 1) and advance hydrodynamic equations 
in time using the Pencil-Code 2 . Because we investigate the 1-D 
case in this work, we cannot apply the 2-D Alternate Direction 
Implicit (ADI) scheme developed in Paper I. We therefore code 
a 1-D version of this scheme that consists of a semi-implicit 



Crank-Nicholson algorithm where nonlinearities are dealt with 
using a Jacobian factorisation (see e.g. Press et al. 1992, §16.6). 



1 Figures are plotted using SciPy, open source scientific tools for 
Python (Jones et al. 2001-2008). 

2 See http://www.nordita.org/software/pencil-code/ and 
Brandenburg & Dobler (2002). 
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Fig. 2. a) Temporal power spectrum for the momentum in the 
(z, 6L>)-plane. b) The resulting spectrum after integrating over 
depth, c) Comparison between normalised momentum profiles 
for n = (0, 2) modes according to the DNS power spectrum 
(solid black line) and to the linear- stability analysis (dotted blue 
line). 
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To determine which modes are present in the DNS once 
it has reached the nonlinear- saturation regime, we perform in 
Fig. 2a a temporal Fourier transform of the momentum field 
pu(z, t) and plot the power spectrum in the (z, 6L>)-plane. With 
this method, acoustic modes are extracted because they emerge 
as "shark fin profiles" about definite eigenfrequencies (see 
Dintrans & Brandenburg 2004). In Fig. 2b, we integrate pu(z, co) 
over depth to obtain the mean spectrum. Several discrete peaks 
corresponding to normal modes appear but the fundamental 
mode close to a>o = 5.439 clearly dominates. 

The linear eigenfunctions can be compared to the mean pro- 
files computed from a zoom about given eigenfrequencies in the 
DNS power spectrum shown in Fig. 2. Figure 2c displays such 
profiles about co^ = 5.439 and 0J2 = 11.06 (solid black lines), 
whereas associated eigenfunctions are overplotted in dotted blue 
lines. The agreement between the linear- stability analysis and 
the DNS is remarkable. 

In summary, Fig. 2 shows that several overtones are present 
in this DNS, even for long times. Because these overtones are 
linearly stable (see Table 1), some underlying energy transfers 
occur between modes. To investigate these intricate couplings in 
more detail, we need to follow the evolution of each mode sepa- 
rately using the projection method developed in the next section. 

3. The adjoint problem 

As said in the Introduction, the radiative diffusivity in our model 
is large close to the surface due to the instability criterion (the so- 
called "transition region" criterion). Eigenvectors strongly dif- 
fer from adiabatic ones, and therefore the quasi-adiabatic ap- 
proximation fails in this case. The dissipative effects must be 
fully taken into account and the oscillations operator is non- 
Hermitian. In other words, the matrix A provided in Eq. (7) is 
not self- adjoint and its eigenvectors T y - are not mutually orthog- 
onal. This implies that they cannot be used to arrange an acoustic 
subspace on which the physical fields of the DNS are projected. 
We thus consider the adjoint problem that has the same spec- 
trum, while its eigenvectors are orthogonal to A-ones because of 
the biorthogonality property (see Appendix A.l). 

3. 1 . The adjoint matrix A 1 " 

We start from the following inner-product 

-1 



(X, Y)= \ X 1 
Jo 



Ydz, 



(8) 



where f denotes the Hermitian conjugate. The adjoint of the op- 
erator A is defined by (e.g. Lowdin 1983) 



(x,ay) = (a^x,y). 



(9) 



The eigenvectors <D; of A^ are normalised to verify (see Appen- 
dix A.l) 



(*t 9 Wj) = 6 iJ9 
where <D; are the eigenvectors of A 1 " such that 



Mj = Wb 



AQ>i = A*Q>i. 



(10) 
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The adjoint-matrix calculation is detailed in Appendix A. 2 
and the resulting A 1 " is provided in Eq. (16). This corresponds to 
the following oscillation equations for the adjoint problem: 
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3.2. The adjoint boundary conditions 

These boundary conditions are detailed in Appendix A. 3. While 
integrating Eq. (A.6) by parts, we obtain surface contributions 
that must vanish to fulfill the adjoint definition (9) (Bohlius et al. 
2007). This implies the following set of adjoint boundary condi- 
tions 



u = 0forz = (0, 1), 

dV (dK dlnp \ 

K — -j- +/ ^o— } — \T' = 0forz = 0, (13) 

dz \ dz dz j 

V =0forz= 1. 



3.3. Results 

We solve the adjoint eigenvalue problem defined by Eqs. (12-13) 
again using the LSB code, that is, we compute both the whole 
spectrum of eigenvalues A* n and their associated eigenvectors Q> n . 
An example is given in Fig. 3 where we plot the real part of 
both the regular (solid black line) and adjoint (dashed blue line) 
eigenfunctions of the fundamental n = mode. 



4. The projection method 

Once the two sets of (regular) *F W and (adjoint) Q> n eigenvector 
are known, the DNS fields are projected at each timestep using 
the following procedure: 

1 . We first arrange the physical fields as 



^DNsfcO 



(T'\ 
u 
R) 



(14) 



DNS 



where u(z, i) is the velocity and T'(z, t), R(z, t) the tempera- 
ture and density perturbations provided by 



T'(z,t) = nz,t)-T (zl R(zj) 



// = p(z, t) 
Po Po(z) 



-1. (15) 
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Fig. 3. Real part of normalised eigenfunctions (u, R, T') for the 
fundamental mode (solid black line) and its adjoint (dotted blue 
line), for the equilibrium setup denoted by a white cross in Fig. 1 . 

2. We decompose these fields onto the regular-Xmzdx eigenvec- 
tors 



E DNS fc t) = % i £ c n (tyP n (z) \ , (17) 



7? = 



where the c n (t) coefficients are completed using adjoint 
eigenvectors 



c n (t) = <O w , S DNS ) = I O^Sdns*- 
Jo 



(18) 



This complex function c n (t) defines the mode amplitude be- 
cause its time evolution is related to the eigenfrequency o) n 
by 



C n (t) = \c n (t)\ e 1M) With (f> n (t) oc OJ n t. 



(19) 



3. We perform Fourier transforms or phase diagrams of these 
c n (t) coefficients to investigate the evolution of modes. 

4.1. Time evolution of the complex-mode amplitudes c n (t) 

The left panels of Fig. 4 show the time evolution of the real 
part of the projection coefficients c n {f) for the first four modes 
n e [0, 3 J. We overplot the curves oc exp(r„ t) derived from the 
linear- stability analysis (see Table 1). As expected, only the fun- 
damental mode amplitude co(t) grows at each timestep. For other 
linearly- stable modes, amplitudes decay except for the n = 2 
one. The amplitude of this mode shows a transient decrease un- 
til time t =* 10 and then begins to increase. This behaviour is 
an estimate of the nonlinear coupling that occurs between this 
overtone and (at least) the fundamental mode. 



We next perform Fourier transforms of these mode ampli- 
tudes (right panel of Fig. 4). Following Eq. (19), the c n (t) am- 
plitudes behave as c n (t) oc exp(k*v) because theoretical eigen- 
frequencies and peaks appearing in power spectra overlap one 
another. 



4.2. Phase diagrams 

In each panel of Fig. 5, the imaginary part of the complex mode 
amplitude c n (t) is plotted as a function of the real part for the 
four modes shown in Fig. 4. To avoid unreadable plots, we stop 
the time evolution at t = 50. The decreasing orbits are plotted 
as black lines, whereas the increasing ones are in blue. In these 
diagrams, the radius of each trajectory is related to the modulus 
of the complex mode amplitude, which implies that a decreasing 
(increasing) radius is associated with a damped (excited) mode. 
Once again, the fundamental mode appears to be continu- 
ously excited. We note that the excitation phenomenon is coher- 
ent in the sense that no discontinuity is observed in the mode 
orbits, that is, the spiral develops continuously. For the second 
overtone, after the linearly-transient decrease, the radius begins 
to increase significantly and this is still the signature of a nonlin- 
ear coupling with the fundamental mode. For the linearly- stable 
n - \ and n = 3 modes, things are more intricate because a 
marginal increase is shown during the last orbits (also seen in 
Fig. 4). Unfortunately, these phase diagrams are not adapted to 
investigate this long-time dynamics because orbits will overlap 
one another. 



4.3. The mode kinetic energy content 

To address the mode couplings responsible for the nonlinear 
limit-cycle stability observed at late times, we could compute 
the energy content of each mode separately. 

In BCM93, the total amount of sound generated by the turbu- 
lent convection is weak because they find that only 0.22% of the 
kinetic energy is stored in acoustic standing waves. Our prob- 
lem however consists of an initially static radiative zone without 
convection and the velocity field that develops is only due to 
acoustic modes, that is, 



F tot _ 



= f>», 



(20) 



72-0 



where E n is the kinetic energy contained in the ^-acoustic mode. 
One advantage of our projection method lies in its ability to com- 
pute this kinetic energy content E n because Eq. (17) implies that 



^dns 



% 



^ C n (t)u n 



(21) 



U=o 



where wdns and u n are the velocity of the DNS and the regular 
eigenvector, respectively. By multiplying by l/2pou, the kinetic 
energy density E^ n (z) is obtained as 
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Fig. 4. Left panel: time evolution of the real part of the complex-mode amplitudes c n (i) for n e [0, 3] modes. The theoretical growth 
or damping curves are overplotted in dashed blue lines. Right panel: their respective Fourier transforms. The modes eigenfrequencies 
are overplotted in dashed blue lines. 
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After integrating over the domain, the total kinetic energy reads 



00 ( °° 1 r^ 

4in=y y« p %[Cn(t)u n ]%[c p (t)u p ]dz 



n=0 \p=0 



(23) 



It is possible to determine the contribution of each acoustic mode 
to the whole amount of kinetic energy because Eq. (20) allows 
us to define E n as 



00 1 M 

p~i 2 J o 



P0%>[c n (f)Unl%>[Cp(f)Up]dz. 



(24) 



Figure 6a displays the long-time evolution of the kinetic en- 
ergy ratio E n /E^ n for n e [0, 6]. After the transient linear 
growth of the fundamental mode, a given fraction of energy 
is progressively transferred to upper overtones and the nonlin- 
ear saturation is achieved above t ^ 150. These nonlinear cou- 
plings mainly involve the n = and n = 2 modes because their 
energy ratios are dominant. However, the other overtones (i.e. 
n g [1,3-6]) also participate in this nonlinear saturation pro- 
cess but their energy ratios remain below 1%. This is compatible 
with the weak peaks around the corresponding frequencies oj n 



shown in the power spectrum in Fig. 2b or with the marginal 
increase of the radius orbits in Fig. 5. It is at first glance not 
surprising because a high-order mode is associated with smaller 
scales and therefore with a high damping rate (Table 1). This im- 
plies that the excitation of these overtones from the fundamental 
mode would require efficient underlying couplings that are not 
observed in this simulation. 



The fact that the n - 2 mode is more excited than the n - 1 
one nevertheless contradicts this explanation because its damp- 
ing rate is more than twice the n — 1 one. In Fig. 6b, we expand 
Fig. 6a for the fundamental mode and this second overtone only. 
The major energy transfer well occurs between these two modes 
because they contain about 98% of the total kinetic energy of 
this DNS. The reason for this favored coupling lies in the period 
ratio existing between these two modes (see Table 1): the funda- 
mental period is Po = 2n/aJo =* 1.155, while the n = 2 one is 
P2 - 0.568 such that the corresponding period ratio is close to 
one half (P2/P0 - 0.491). The n = 2 mode therefore receives in 
a preferential way some energy from the n = mode every two 
periods and that corresponds, from a dynamical point of view, to 
a classical 2:1 resonance. Such a resonance is usual in celestial 
mechanics with, e.g., Jupiter's moons Io (P = 1.769d), Europa 
(P = 3.55 Id) and Ganymede (P = 7.154d) and it is well known 
that it helps to stabilise orbits. In our case, this stabilisation takes 
the form of a nonlinear saturation: the linear growth of the fun- 
damental mode is balanced by the pumping of energy from the 
linearly- stable second overtone behaving as an energy sink. The 
final amplitudes give about 87% of the kinetic energy in the fun- 
damental mode and 1 1 % in the second overtone, the remaining 
2 percent being in higher overtones as displayed in Fig. 6a. 
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Fig. 5. Phase diagrams for n e [0, 3] modes. The abscissa (ordinate) corresponds to the real (imaginary) part of the complex mode 
amplitude c n (t). The increasing orbits are plotted in blue, whereas the decreasing ones are in black. 



5. The Hertzsprung progression 

The nonlinear saturation in our excitation model therefore re- 
sults from a 2:1 resonance between the fundamental mode 
and the second overtone. As shown in the Introduction, this 
mechanism has also been proposed to explain the secondary 
bump observed in the luminosity variations of Bump Cepheids 
(SS76). An interesting correlation between the value of the ra- 
tio P2/P0 and the bump position was also emphasised leading 
to the so-called "Hertzsprung progression" (Simon & Lee 1981; 
Kovacs & Buchler 1989; Buchler et al. 1990). Nevertheless this 
correlation with P2/P0 is not observed as the period P2 deduced 
from luminosity curves is necessarily equal to Po/2. Indeed, the 
second overtone period changes due to nonlinear couplings and 
finally reaches the value P2 = Po/2 once the nonlinear saturation 
is achieved (i.e. the 2:1 resonance). As a consequence, only lin- 
ear eigenperiods are relevant when studying the influence of the 
P2/P0 ratio on the bump position. Furthermore, Buchler et al. 
(1990) have shown that this correlation between the phase and 
the P2/P0 ratio is not reproduced in their models when plotting 
the phase as a function of Pq only. As our hydrodynamical model 
deals with dimensionless quantities, studying the phase varia- 
tions according to P2/P0 is also more consistent. 

To determine whether our model is able to reproduce this 
Hertzsprung progression, we first overplot in Fig. 7 the P2/P0 
level lines with the instability strip discovered in the linear- 
stability analysis (Section 2.1). The solid white level lines corre- 
spond to P2/P0 > 1/2, the dashed lines to P2/P0 < 1/2 and the 
red line is for the critical level P2/P0 - 1/2. The instability strip 
is clearly split in two separate parts and the most unstable modes 
belong to the P2/P0 < 1/2 region. If we now perform the DNS 
corresponding to those different setups, we can check if the lu- 
minosity bump position is really correlated with the period ratio 



value. Because we are dealing with the 1-D case, the luminosity 
reduces to the radiative flux at the top of the domain, that is, 



L(t) = Ftop(r), 



(25) 

and the luminosity curves are transformed into emerging radia- 
tive flux curves. Because the bottom flux and the top tempera- 
ture are imposed in the DNS, these flux variations F top (t) happen 
about Fbot- 

We then run a large number of DNS (up to 400) in the in- 
stability strip to cover a significant range in the P2/P0 ratio and 
cross the 1 /2-line; all these DNS belong to the dashed box shown 
in Fig. 7. In each case, the simulation is performed above the 
nonlinear saturation fend = 500). Because we know that the fun- 
damental mode and the second overtone enter in the time evolu- 
tion of the emerging flux F top (t), we fit the obtained flux varia- 
tions using the following decomposition: 

fit) = F hot + A cos(oj t + 0o) + A 2 cos(2oj t + 2 ), (26) 

where the free parameters [Ao, A2, <po, 02 ] are computed us- 
ing the Levenberg-Marquardt algorithm that is a nonlinear least- 
square method (e.g. Press et al. 1992, §15.5). The phase shift 
between the two signals is defined by 



A0 = 02 - 200, 



(27) 



and allows first to phase the flux curves and second, to know 
where the secondary bump is. If A0 is less (greater) than n, this 
secondary bump will be located in the descending (ascending) 
branch after (before) the flux maximum. 

The whole DNS are summarised in Fig. 8 where isocontours 
of A0 are plotted in the same (7b U mp, &0 instability strip plane 
used in Fig. 7. The red regions are associated with A0 > n, the 
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Fig. 6. a) Kinetic energy ratios for n e [0, 6] in a logarithmic y- 
scale. b) Zoom in the n = (solid black line) and n = 2 (dashed 
green line) mode energy contents. 
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Fig. 7. P2/P0 level lines overplotted on the instability strip. 
Dashed white lines denote ratio P2/P0 < 1/2, while solid white 
lines correspond to P2/P0 > 1/2. The critical level P2/P0 = 1/2 
is emphasised as a red line, and the dashed box delimits the 
region where a large number of DNS are performed to study 
Hertzsprung's progression (see Figs. 8-9). 

blue ones with A0 < n and the white ones with A0 = n. By over- 
plotting the level lines of the period ratio P2IP0, the following 
correspondence appears: 



A0 < 71 (bump after the max.) o P2/P0 > 1/2, 
A0 > 71 (bump before the max.) o P2/P0 < 1/2, 




it 






2.25 2.20 2.15 2.10 2.05 2.00 1.95 1.90 



^bump 



(28) 



Fig. 8. Isocontours of A0 in the (7b um p, ffi) plane (note that we 
rather plot Acp-nto highlight the critical line A0 = n). This fig- 
ure corresponds to the dashed box in Fig. 7 and the black circles 
display eleven selected simulations showed in Fig. 9. 



which corresponds to the observed correlation (SS76). 

To highlight this result, examples of light curves shaping a 
Hertzsprung progression are provided in Fig. 9. These eleven 
curves correspond to the black circles distributed across the 
P2/P0 = 1/2 level line in Fig. 8. The correspondences (28) ap- 
pear clearly: the bump is located in the descending branch of the 
light curve if P 2 /Po > 1/2, while values P 2 /Po < 1/2 lead to a 
bump in the ascending branch. 

Observations also show that the amplitude of the surface ve- 
locity presents a double-peaked behavior along the Hertzsprung 
progression as its center corresponds to a minimum value for the 
velocity (Bono et al. 2000). In our simulations, no clear correla- 
tion between the amplitude ratio A2/A0 and the bump position 
appears and we guess that it may be related to the two following 
restrictions: 

- Rigid walls for the velocity have been assumed. As we are 
interested in the physics involved near the top of the box, the 
ratio A2/A0 is not really relevant there compared to, say, the 
phase difference between the two modes. 

- The sound speed contrast across the box is relatively weak. 

As the frequency step is given by Ay ~ [2 j L dr/c s ) , there 
are not enough variations to obtain significant changes in this 
ratio. 

However our model is able to qualitatively reproduce the 
Hertzsprung progression, as the position of the secondary bump 
with regard to the main bump evolves in the expected way with 
the P2/P0 ratio. We therefore conclude that the 2:1 resonance 
supplying the nonlinear saturation in our model results in an 
"observational" bump in the light curve of which the position 
is linked to the period ratio ^2/^0- 

6. Conclusion 

In this paper, we investigate the nonlinear saturation of the 
acoustic modes excited by the ^-mechanism in direct numerical 
simulations (DNS). This study follows on from our former work 
where we began our modelling of the K-mechanism in Cepheids 
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Fig. 9. Rephased time evolutions of the surface radiative flux 
about time t = 500 for the eleven selected simulations plotted 
as black circles in Fig. 8. The ratio of eigenperiods P2/P0 is in- 
dicated in each case. 



tal mode and the second overtone. This last mode, which rep- 
resents about 10% of the total kinetic energy, is involved in the 
nonlinear saturation of the instability through a 2:1 resonance 
with the fundamental mode. It means that the unstable funda- 
mental mode gives energy to the stable second overtone that 
leads to the limit-cycle stability, as displayed in Fig. 6. 

This 2:1 resonance is striking because the same mecha- 
nism is probably responsible for the observed bump in the lu- 
minosity curves of Bump Cepheids, leading to the well-known 
"Hertzsprung's progression". This progression links the bump 
position to the period ratio P2/P0 existing between the funda- 
mental mode and the second overtone (SS76). We perform a 
large number of DNS covering a significant range in this ratio 
and extract the resulting luminosity curves. The 2:1 resonance 
does lead to an Hertzsprung progression as the bump arises in 
the ascending branch for P2/P0 < 1/2, while it is located in the 
descending branch for P2/P0 > 1/2. 

The study of the A:-mechanism in a purely-radiative layer 
is now completed with this work. The physical criteria for 
the instability as well as the nature of its nonlinear saturation 
are addressed. Future works will be devoted to the influence 
of convection on the mode stability. Indeed, a coupling be- 
tween the convection and pulsations is suspected to play a ma- 
jor role in the disappearance of unstable modes close to the 
red edge of Cepheids' instability strip. Time-dependent theo- 
ries of convection have difficulty in explaining the red-edge lo- 
cation since they rely on many unconstrained parameters (e.g. 
Wuchterl & Feuchtinger 1998; Yecko et al. 1998). Because DNS 
fully account for crucial nonlinearities, they are appropriate to 
properly investigate this dynamical coupling between the acous- 
tic modes and convection. We thus plan to perform 2-D DNS of 
the /c-mechanism in the presence of convection. The interesting 
cases will of course correspond to a (local) convective-turnover 
timescale of the same order as the period of the fundamental 
mode. 
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Interuniversitaire de Calcul de Toulouse" (http://www.calmip.cict.fr/) 
and Grid' 5000, which is an initiative from the French Ministry of Research 
through the ACI GRID incentive action, INRIA, CNRS and RENATER and other 
contributing partners (see https : //www . grid5000 . f r). 



by looking the propagation of acoustic modes in a layer of gas 
in which the ionisation is modelled by a hollow in the radiative 
conductivity profile (Gastine & Dintrans 2007, Paper I) 

To catch the evolution of linearly-unstable modes, we ap- 
ply an original method consisting of the projection of the DNS 
fields onto the eigenmodes that are solutions to the linear os- 
cillation equations (Dintrans & Brandenburg 2004). Because the 
oscillation operator is non-Hermitian, the regular eigenvectors 
are not mutually orthogonal and we consider the adjoint prob- 
lem that has the same spectrum while its eigenvectors are or- 
thogonal to the regular ones. Thanks to these two sets of eigen- 
vectors, we compute the projection coefficients for each DNS 
snapshot and follow the temporal evolution of acoustic modes 
separately. Corresponding phase diagrams show that the orbit 
radius of the linearly -unstable fundamental mode is continu- 
ously growing with time. On the contrary, orbits associated with 
the linearly-stable overtones exhibit decreasing radii, except the 
n = 2 one that is increasing at late times. 

We derive the kinetic energy contributions of each mode 
propagating in the DNS. More than 98% of the total kinetic en- 
ergy is contained in two modes corresponding to the fundamen- 



Appendix A: The adjoint formalism 

A. 1. The biorthogonality relation 

We consider a continuous linear operator H with eigenvalues X\ 
and eigenvectors ¥,- linked by 



HVi = Wi. 



(A.l) 



The operator H^ is its adjoint with eigenvalues /// and eigenvec- 
tors Oj as 



Using the selected inner product (8), we have 

<Yy,i^*,>=/i^,**>, 

that implies 



(A.2) 



(A.3) 



(HV j9 *i) = IMWj.Oi), (A.4) 

by using the adjoint definition (9). We thus obtain 
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(A.5) 



While letting aside in a first step the integrated terms, the three 
terms (a), (b) and (c) imply that 



(A* j -fi i )?¥j,<S> i ) = hence <¥,-,**> = % 

This property is the biorthogonality satisfied by both the regular (a\ + (h\ + (r\ - f T' ^° in 2 -2D In 

and adjoint eigenfunctions (e.g. Lowdin 1983). 7 J 1 p I 



A2. The derivation of the adjoint matrix A 1 " 

The relation (9) that links an operator with its adjoint takes the 
following form in our problem 



-D 2 lnpo + (D lnpo) 2 T£dz 



= (H^T^T[). 



(A. 10) 



Applying Eq. (A. 10), we finally get 



f (T' 2 u 2 R 2 )*A 
= f (T[ Ml tf^A 1 






dz 



(A.6) 



lj = 2^ f/) 2 _ 2D In poD-D 2 In p + (£> lnp ) 2 l (A. 11) 
Po l J 



^2 ) 



dz, 



• Calculation of A 



... T . Al comes from the following inner-product 

where iJ^u^Ri) 1 is an adjoint eigenvector and (7; wi 7?i) J is a Zi 

regular one. To determine the adjoint operator A 1 ", we calculate x 

each element successively from the regular matrix in Eq. (7). (j> A\iU\) = \ T'*\ ot - (y - \)Tc\D 

We provide below an example of such a calculation for the first " ~ J ^ [ Kq 

column of the adjoint matrix given in Eq. (16). 



u\dz. (A.12) 



Integrating by parts implies that 



• Calculation of A^ 



Ajj comes from the following inner-product 

/ = {T' 2 ,A n T[) 

= \ T? — \K D 2 + 2DK D + D 2 K ] T[dz. 
Jo Po L J 

This integral is split in 3 parts as 



Jo po dz 2 Jo Po dz dz 



(T'AnUi) 



(r-i)r D + (2-r) 



-(r-i)K*r ^] . 



^0 



T' v ui 



(A.13) 



(A.7) Using Eq. (5), the integrated part cancels and we obtain 

v ^bot 



A^=(y-l)r Z) + (2- r )-^. 

^o 



(A. 14) 



• Calculation of A^ 



(i) 



(2) 



Because A13 = in Eq. (7), we have 

<r£,Ai 3 i?i> = 0, 



(A. 8) an d therefore 



Jo 



Po dz 



.2 M c 



UD 



AJ 1= 0. 



(A.15) 



(A.16) 



(1) and (2) are integrated by parts, leading to 

(1) = 



T^K dT[ 



Po <fe 
Jo dz 2 



dz\ po I l 



A.3. The boundary conditions satisfied by A 1 " 

To determine the boundary conditions for the adjoint problem, 
we follow a method used in Bohlius et al. (2007). Eq. (A.9) con- 
tains the boundary conditions on temperature because the defi- 
nition of the adjoint operator (9) imposes that every integrated 
term vanishes as 



1 d 2 ( T£Kp 

Po 



T\dz, 



(1) + (2) 



</.) 



(A.9) 



T;*K q dT[ 



Po dz 



+ 2 



(2) 



T 2 dK I 1 „ f 1 J (T?dK \ , , 

-2-— ^r -2 _ -L_£ r<fe. 

Po <fe Jo Jo dz\po dz j 



dz\ po ) l 



T JL^r 

Po dz 1 



0. 



(A. 17) 



(c) 



Applying Eq. (5), this equation implies the two following bound- 
ary conditions on temperature: 
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- At the top z = 1: 



t;*k dT\ 

-^-^—±=0^T'=0. 

Po dz 



- At the bottom z = 0: 



'T£_dKo d (T£Ko 



po dz dz \ po 
which simplifies to 



T[ = 0, 



dT 2 idK dlnp \ 



(A.18) 



(A.19) 



(A.20) 
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